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ABSTRACT 


The  problem  of  rapidly  and  accurately  determining  the  orbital  plane  of 
an  artificial  earth  satellite  from  a short  series  (<^  10m)  of  observations  is 
posed  and  solved.  The  observations  can  be  either  angles  only  data  or  angles 
and  angular  rates.  The  accuracy  Is  generally  5°.  In  addition,  another  look 
is  given  to  the  angular  velocity  distance  estimation  technique. 
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I.  MOTIVATION 

During  searches  of  the  equatorial  belt  one  will  find  known  objects 
and  unknown  ones.  For  the  known  ones  immediate  identification  is  desirable 
so  that  telescope  time  can  be  freed  for  other  purposes.  Currently  this  is 
done  by  position  only.  The  use  of  a spin  period,  angular  velocity,  or  an  I 

orbital  plane  direction  would  help.  Hence,  it  was  asked  whether  or  not 
an  artificial  satellite's  orbital  plane  could  be  accurately  determined  in  a 
relatively  short  (5  - 10m)  time. 

If  an  object  turned  out  to  be  an  unknown,  then  finding  it  again  in  the 
immediate  future,  so  that  the  telescope  could  be  time  shared,  was  a problem. 

For  low  eccentricity,  low  inclination  objects  the  current  dead  reckoning 
tools  are  adequate.  For  high  eccentricity  or  high  inclination  objects  these 
techniques  decay  in  predictive  ability  as  the  square  of  the  elapsed  time. 

If  one  knew  the  orbital  plane  then  one  could  constrain  the  predictions  of 

right  ascension  and  declination  to  satisfy  the  orbital  plane  equation  * 

[cf.  Eq.  (4b)].  Once  again  a rapid,  accurate  determination  of  the  orbital 
plane  is  desirable. 

Finally,  there  is  the  problem  of  recovering  unknown  satellites  the  next 
night.  Obviously  an  along  the  orbit  search  needs  the  orbital  plane  to  proceed. 

Also,  discriminating  against  several  possible  initial  orbits  would  be  helpful. 

Hence,  the  need  to  determine  the  orbital  plane. 

Below  we  examine  the  physics  and  geometry  of  the  problem  and  test  a 
large  number  (34)  of  different  artificial  satellite  orbits  based  on  the 
assumption  they  were  discovered  during  equatorial  search.  The  results  are 
generally  good  to  5°. 
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II.  THE  CONSTRAINT  r ' L - 0 

A unit  vector  in  the  direction  of  the  angular  momentum  vector  is 
Riven  by 

L - (sinftsini , -cosftsim  , cosl),  (1) 

where  ft  is  the  longitude  of  the  ascending  node  on  the  equator  and  \ is 
the  inclination  of  the  orbital  plane  to  the  equator.  fte[0,360),  ie[0,180), 

and  the  zero  of  ft  is  the  vernal  equinox.  The  magnitude  of  the  angular 

2 2 1/2 

momentum  vector  per  unit  mass  is  given  by  na  (1  - e ) where  n is  the 
mean  motion,  a is  the  semi-major  axis,  and  e is  the  eccentricity.  Since 

L - r x v/|r  x v| , (2) 

where  r ■ (x,y,z)  - r(cos6cosa,  cos6sina,  sin6)  is  the  geocentric  location 
and  v - dr/dt  is  the  geocentric  velocity,  it  follows  that 

jr  * L ■ 0 and  v ' L = 0.  (3) 

These  relationships  are  true  at  all  instants  of  time  when  a,  e,  l,  and  ft 
are  interpretated  as  the  elements  of  the  osculating  ellipse.  Hence,  no 
assumptions  about  the  forces  acting  on  the  object  (an  artificial  satellite 
bound  to  the  earth)  need  be  made  except  that  di/dt,  dft/dt , da/dt,  and 
de/dt  multiplied  by  the  duration  of  the  observing  span  be  very  small.  It 
so  happens  that  for  both  practical  reasons  concerned  with  the  mechanics  of 
equatorial  search  and  analytical  ones  discussed  below  the  observing  span 
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will  be  on  the  order  of  a few  minutes.  Hence,  we  can  take  i,  ft,  a,  and  e 
to  be  constants. 


\ 
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The  first  part  of  Eq.  (3)  can  be  written  as 


xsinftsini  -ycosftsint  + zcosl  = 0,  (4a)  | 

or,  in  spherical  coordinates, 

r [cos6tanxsin(a  - ft)  - sin6]  - 0.  (4b) 

Note  if  l = 90°  the  right  ascension  is  always  equal  to  ft.  We  also  note  that 
the  geocentric  distance  plays  no  real  role  in  Eq.  (4b).  In  addition,  two 
observations  of  geocentric  position  (e.g.,  a and  6)  suffice  to  determine 
1 and  ft,  viz, 

j 

sinot.  tan5„  - sina„tan6,  „ „ 

tanft — ^ ± -ji*  <5a> 

cosa^tan^  ~ cosot2tan^ 

tarn  = tan61csc(a1  - ft)  - tan62csc(a2  - ft).  (5b) 

The  problem  is  that  we  measure  topocentric  position  A,  A [the  topocentric 
location  is  given  by  R = R(cosAcosA,cosAsinA,  sinA) ] and  not  the  geocentric 
values . 

Let  the  observer's  geocentric  location  be  given  by  s = s(cosif>’cosT, 
cosifi'sinT,  sintfi ' ) where  s is  his  geocentric  distance,  <J>'  is  his  geocentric 

i | 
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j 
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latitude,  and  r the  local  sidereal  time.  Then  the  relationship  between 
geocentric  and  topocentrlc  locations  Is 

r - R + s . (6) 

In  terms  of  the  measured  quantities  Eq.  (4b)  becomes 

R[cosAtanlsin(A  - ft)  - slnA] 

+ s[cos<J>’tanisln(T  - ft)  - sin$']  ■ 0.  (7) 

We  see  that  the  distance  does  not  now  drop  out.  Moreover,  if  it  happens 
that  even  within  the  5 - 10m  observing  span  R changes  appreciably,  Eq.  (7) 
will  not  accurately  reflect  the  physics  of  the  problem.  Hence,  we  allow 

for  a constant  variation  of  R,  V ■ dR/dt.  Therefore,  r ’ L • 0 becomes 

f(\,  ft,  R,  V,  t)  - (R  + Vt){cosA(t)tan\sin[A(t)  -0]  - 
sinA(t)}  + stcos^’tanlsinfTCt)  - ft]  - sin$'}  ■ 0.  (8) 

The  full  time  dependence  is  now  explicitly  indicated  in  Eq.  (8). 

Presumably  we  have  N > 3 measurements  of  A,  A,  and  t,  say 

{ Aj , Aj , tj},  i ■ 1,  2,  3,  ...,  N.  From  these  and  Eq.  (8)  we  wish  to 
determine  i and  ft  and  te  apparently  get  R and  V free.  Letting 
fj(l,  ft,  R,  V)  - f(t,  ft,  R,  V,  tj)  we  form 
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[i,  si,  r,  v)  - l w f 2(i,  si,  r,  v) 

j-i  J J 


and  seek  the  minimum  of  F with  respect  to  \,  SI,  R,  and  V.  The 
quantities  w^  are  the  weights  of  the  jth  value  of  f.  We  argue  below 
for  Wj  - 1/N  Vjell.N] . The  problem  is  now  in  the  form  of  a non-linear 
least  squares  estimation  problem.  There  exists  a standard  technique 
for  solving  such  problems  (see  the  Appendix).  Unfortunately,  it  doesn't 
work  in  this  case.  Instead  the  method  of  steepest  descent  has  been 
successfully  used  to  obtain  values  of  \,  SI,  R,  and  V (see  the  Appendix). 
The  accuracy  of  the  results  obtained  by  using  Eq.  (9)  is  discussed  in 


The  weights,  w^ , are  given  by 


i/Wj  - ofj/aAj)2^  + ofj/aAj)2^  + ofj/atj)2^  , qo) 


where  w is  the  weight  of  A , etc.  Relative  to  the  weights  of  topocentric 

J j 

right  ascensioi  *nd  declination,  the  weight  of  the  time  is  effectively 
infinite.  Also,  since  A^  = 0°  Vje[l,N],  the  weights  of  right  ascension 
and  declination  are  effectively  equal.  Finally,  as  t^  - t^  ~ 10™,  to  first 
order  all  the  3fj/3Aj  and  3fj/9Aj  are  the  same,  whence  to  first  order,  all 
of  the  weights  of  the  same.  After  reading  the  Appendix  it  will  become 
clear  that  a dependence  of  the  weights  on  l,  SI,  R,  or  V is  immensely 
complicating  analytically.  Hence,  even  if  it  weren't  true  that  the 


kL 


weights  are  all  equal  (to  first  order),  1 would've  taken  them  to  be  equal 
out  of  computational  necessity.  This  discussion  also  makes  it  clear  that 
determining  11,  and  R must  rest  on  second  order  effects  so  that  high 
accuracy  (e.g.,  better  than  1°  In  L)  is  not  to  be  expected.  Clearly, 
determining  V rests  on  third  order  effects. 
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III.  THE  CONSTRAINT  v - L - 0 

Unless  r x v vanishes  (e.g.,  e - o)  the  constraint  v ' L - 0 is 
independent  of  the  constraint  r ’ L - 0.  Assuming  that  the  acceleration 
in  the  topocentric  range  is  negligible,  we  can  write  this  as 

df(l,  ft,  R,  V,  t)/dt  - V{cosA(t)tamsin[A(t)  - ft]  - sinA(t)} 

- (R  + Vt)A  (t){sinA(t)  tan\sin[A(t)  - ft]  + cosA(t)} 

+ (R  + Vt)  A(t)cosA(t)  tanlcos(A(t)  - ft] 

+ sx(t)cos4>'tanicos[T(t)  -ft]  *0.  (11) 

The  function  to  be  minimized  is  now 

N . 

f(i.  ft,  r.  V)  > l w [f,(\.  n,  R,  V)  + cU\,  n,  r,  v)].  (12) 

J-l  J J 1 

The  weights  for  f and  df/dt  are  taken  to  be  same  because  the  f weights 
are  all  equal  to  each  other  to  first  order.  We  note  that  no  new  unknowns 
are  introduced  in  Eq . (12)  over  those  already  in  Eq.  (9). 

Solving  Eq . (12)  is  no  more  difficult  than  solving  Eq . (9)  and  it 
has  also  been  performed  by  the  method  of  steepest  descent.  If  it  turns 
out  that  no  information  can  be  obtained  about  V even  when  it's  relatively 
large  we  merely  set  V = 0 in  Eqs.  (8  and  11). 
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IV.  NUMERICAL  TESTS 


The  time  spacing  for  the  first  series  of  tests  was  tj+j  ~ tj  = 30 
Vje[l,N] . This  series,  using  six  widely  different  artificial  satellite 


orbits,  included  V,  used  the  method  of  steepest  descent,  and  showed  that 
for  a =5m  observing  span  the  direction  of  L could  be  gotten  within  a few 
(5)  degrees.  Eliminating  V from  the  equations  made  no  difference  (it  was 
always  O.OOOkm/sec  anyhow)  even  though  |v|  was  as  large  as  5km/sec  in  some 
cases.  I also  rewrote  Eqs.  (8  and  9)  using  f/R  (see  the  Appendix).  For 
these  same  satellites  and  the  same  observations  the  results  were  never 
better  (with  or  without  V)  and  sometimes  poorer.  Thus,  it  seemed 
important  to  test  a wider  variety  of  satellite  orbits,  the  effect  of  the 
angular  rates,  the  effect  of  the  quality  of  the  data,  and  the  effect  of 
the  time  span  of  the  data. 

A total  of  34  different  satellite  orbits  have  been  used  to  generate 
topocentric  position  and  topocentric  angular  velocity.  Each  has  been 
analyzed  (without  V)  in  the  following  ways;  (a)  t^+^  ~ tj  = 30s  V j e f 1 , N ] , 
accuracy  in  position  = \/2'' , N = 11  (e.g.,  5m)  , no  angular  velocity,  (b) 
same  as  case  (a)  except  N = 21  (e.g.,  10m) , (c)  same  as  case  (a)  but  with  a 
positional  accuracy  of  9"  (SSCSAO  quality),  (d)  same  as  case  (c)  except 
N = 21,  (e)  same  as  case  (a)  except  including  angular  velocities  accurate  to 
O'.'l/sec,  and  (f)  same  as  case  (e)  except  N = 21. 

If  w and  I are  our  approximations  for  and  i,  then  let 


£.  = (sinwsinl,  -coswsinl  ,cosl)  , 
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If  8.  and  L were  parallel  then  <f>  would  be  zero.  As  there  is  an  ambiguity  In 


determining  ft  from  F.q.(5a),  we  may  find  a set  of  values  for  I and  w which 


has  2.  anti-parallel  to  L.  This  corresponds  to  (J>  ■ 180°  and  the  Incorrect 
heliticity.  Lastly  we  note  If  i and  I are  both  small,  may  be  small  but 
|w  - ft|  not  necessarily  small. 

Within  a few  tenths  of  a degree  in  <p  cases  (a),  (b) , (c) , and  (d)  all 
produce  the  same  results.  Table  1 lists  the  case  (c)  and  case  (e)  results. 
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The  satellites  are  grouped  by  increasing  mean  motion  and  within  each  mean 


motion  group  by  increasing  eccentricity.  In  general  adding  the  angular 
velocity  information  produces  better  results  (as  would  be  expected).  In  a 
few  instances  the  reverse  is  true.  Moreover,  one  could  presumably  find  a 
function  of  the  satellite's  orbital  elements  and  the  observer's  location 
which  would  predict  the  values  of  <f>  in  Table  1.  Whereas  this  would  explain 
Table  1,  it  is  obviously  of  no  practical  use. 

If  we  concentrate  on  those  satellites  with  a mean  motion  of  less  than 
2.5  rev/day  the  mean  value  of  <f>  is  5!8.  The  standard  deviation  about  the 
mean  is  5 ? 3 . It,  therefore,  appears  that  the  usefulness  of  these  results 
is  fairly  limited.  Further  experimentation  for  dead  reckoning  purposes 
has  been  temporarily  suspended.  Finally,  by  examining  the  results  for  R, 
it  became  clear  that  in  the  huge  majority  of  the  cases  what  came  out  was 
what  went  in.  Hence,  a more  detailed  examination  of  the  angular  velocity 
distance  estimation  technique  seemed  in  order. 
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V.  THE  DISTANCE  ESTIMATION  TECHNIQUE* 

As  this  is  the  largest  body  of  data  I have  yet  had  with  which  to 

test  the  angular  velocity  distance  estimation  technique,  I performed 

such  an  analysis.  Figure  1 is  a plot  of  the  estimated  distance  against 

the  true  distance  (both  are  geocentric).  Low  eccentricity  orbits  are 

plotted  as  dots,  high  eccentricity  orbits  as  plus  signs.  The  two  lines 

represent  r - r and  r - 0.83r.  The  latter  includes  the  correction 
est  est 

of  the  unknown  eccentricity  assuming  a uniform  distribution  of  eccentricities 
in  the  artificial  satellite  population.  Once  beyond  = 3 earth  radii  the 
technique  works  quite  well.  Of  course,  for  all  the  small  eccentricity 
objects  regt  is  systematically  too  low  by  1 - 0.83  - 17%.  To  remedy 
this,  since  about  75%  of  the  total  deep  space  artificial  satellite 
population  is  low  eccentricity,  we  can  split  the  difference  (e.g.,  use 
0.91  instead  of  0.83)  or  use  a two  point  eccentricity  distribution.  The 
latter,  if  75%  are  near  e » 0 and  25%  near  e - 0.7  yields  0.94.  In  the 
future  we'll  split  the  difference. 


*See  SVIIA  of  L.  G.  Taff,  "Astrometry  in  Small  Fields,"  Technical  Note 
1977-2,  Lincoln  Laboratory,  M.I.T.  (14  June  1977),  DDC  AD-A043568/5 . 
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Fig.l.  Plot  of  estimated  geocentric  distance  (rest)  vs  the  actual 
geocentric  distance  (r).  Plus  signs  denote  high  eccentricity 
objects.  Dots  denote  low  eccentricity  objects. 


APPENDIX 


In  this  section  I want  to  review  the  standard  technique  for  solving  a 
non-linear  estimation  problem,  review  the  method  of  steepest  descent,  and 
examine  the  application  of  these  techniques  to  Eqs.  (9  and  12).  Suppose  we 
have  some  function  g of  the  unknown  parameters  x and  we  wish  to  minimize 


G = 1 w g^(x). 

J-l  J J 


(Al) 


The  standard  technique  to  find  the  value  of  x,  say  x^,  that  minimizes  G 


is  to  expand  g (x)  about  some  guess  for  x , say  x . Then,  to  first  order, 
j — — m — o 


G VW  + Vj-^x  -2^’  (*  - V]2- 


(A2) 


G,  as  given  by  Eq.  (A2)  is  now  a linear  function  of  x - and  we  form 

the  normal  equations  as  in  the  ordinary  linear  case  (e.g.,  demand 

V G = 0) . After  solving  the  normal  equations  for  x - x our  new  guess 
X o 

for  x is  x = (x  - x ) + x . If  necessary  this  iteration  is  repeated 

— m — n — — o — o 

until  |x  - x^J  is  sufficiently  small  and  |VxG|  is  sufficiently  small  too. 

This  procedure  won't  work  for  Eqs.  (8  and  9)  because  of  the  following: 

It  is  logical  to  use  the  geocentric  relationships,  Eqs.  (5),  to  obtain 

guesses  for  i and  ft.  We  can  guess  R by  the  angular  velocity  technique 

and  V * 0 should  do.  The  matrix  of  the  normal  equations,  along  the  main 

diagonal,  will  then  have  terms  proportional  to  N [cosAtanisin(Aj  - ft) 

2 2 

- sinAj ] and  (t^  - t^)  times  this.  Clearly  each  of  these  is  very  small, 
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the  second  being  10  ^of  the  first.  Hence,  inversion  of  the  matrix  is 

numerically  impossible.  One  might  think  that  this  could  be  avoided  by 

using  f/R  instead  of  f in  Eq.  (9).  The  present  difficulty  would  be 

removed  but  a new  one  introduced.  As  s/R  = 1/6  for  the  satellites  of 

interest,  the  matrix  of  the  normal  equations  now  contains  terms  of  the 

order  of  N(s/R  ) and  N(t^  - t^)  (s/R  ) on  the  main  diagonal.  These  are 

-2  -7 

of  the  order  of  10  and  10  , respectively.  Hence,  numerical  inversion 
is  again  impossible. 

Having  exhausted  analytically  rearranging  Eqs.  (8  and  9)  we  turn  to  a 
completely  different  technique.  Let  us  return  to  Eqs.  (A1  and  A2)  but 
continue  the  expansion  of  G(x  + AX)  to  second  order  terms.  If  T (x)  is  the 
Hessian  matrix  of  G then 

G(x  + AX)  - G(x)  + VxG(x)  • AX  + (A2/2)  X ‘ T(x)  * X.  (A3) 

As  we  are  looking  for  a minimum  of  G we  take  X to  be  in  the  direction  of 

the  maximum  rate  of  decrease  of  G,  i.e., 

X=-VxG(x).  (A4) 

With  this  value  for  X we  can  find  the  value  of  A to  use  by  insisting  that 
G[x  - AVxG(x)]  be  a minimum  with  respect  to  A.  The  result  is 

A - |VxG(x)|2/[VxG(x)  ' r(x)  • VxG(x)l.  (A5) 
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Thus,  starting  from  a guess  x - x , for  x , we  compute  our  new  guess  from 

— — o —tn 


x - X - XV  G(x) I 
-n  -o  x — ' x ■ x 


with  X evaluated  via  Eq.  (A5)  at  x • x^,-  If  G(x)  has  a unique  minimum 
(e.g.,  x^)  In  any  closed  region  and  the  metric  defined  by  T(x)  has  a 
positive  upper  bound  In  this  region  then  this  iteration  does  converge  to 
x^.  For  most  practical  problems  this  technique  is  difficult  to  apply 
because  one  needs  accurate  values  for  T(x). 

For  the  problem  posed  by  Eqs.  (8  and  9),  or  more  generally  Eqs. 

(11  and  12),  all  of  the  second  derivatives  can  be  computed.  In  particular, 
for  Eqs . (11  and  12) , 


3F/3x  - 2 ^WjIfjOfj/30  + fjOfj/31)] 


etc. , and 


- 2 l WjIOfj/31)2  + fJ(32fj/3i2)  + (3fj/3i)2  + f ^ (32f  j/312)  ] , (A8a) 


3^/3130  - 2 l w [3f  ,/3i)(3f  ./3fi)  + f(32f  /3\3fl)  + (3f  /3i ) (3f  /3Q) 
J J J J J J J 


+ f «‘f  /3i3fl)], 


(A8b) 
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The  needed  derivatives  of  and  df^/dt  are  given  by 


9fj/9l  - sec  x [ (R  + Vt j)cosAjSin(Aj  - ft)  + scos^'sin(T ^ - ft) ] , (A9a) 


9 f ^ / 3ft  - -tanx[(R  + Y't ^ )cosA^cos(A^  - ft)  + scosi|>'cos  (T  j - ft)],  (A9b) 


9f^/3R  - cosAjtan\sin(Aj  - ft)  -sinA^ , (A9c) 


afj/av  » tjDfj/aR,  (A9d) 


3fj/3i  * sec* x [VcosA ^sin(A^  - ft)  - (R  + Vt ^ )A^sinAjSin(A^  - ft) 


+ (R  + Vt^ )AjCosAjCOs(A^  - ft)  + sTjCos^'cos(Tj  - ft)),  (AlOa) 


3f  /3ft  “ -tanx [VcosAjCOs (A^  - ft)  - (R  + Vt ^ )A^ s inA ^ cos (A^  - ft) 


(R  + Vt  j)A^cosA^sin(Aj  - ft)  - st^  cost])'  sin  (t^  - ft)],  (AlOb) 


3f }/3R  * -Aj [sinAjtanxsin(A^  - ft)  + cosA^ ] 


+ AjCOsAj tanxcos (Aj  - ft). 


3f?/av  - 3f  /9R  + t^fj/aR, 


(AlOc) 


(AlOd) 
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32f j / 3 X 2 - 21811x3^/31, 

(Alla) 

32fj/3i3ft  - sec2xcotx 3f ^ /3fl, 

(Allb) 

32fj/3\3R  ■ sec2\cosAjSln(Aj  - ft). 

(Allc) 

32ff/3\3V  - tJ32fJ/3l3R, 

(Alld) 

"fj/3ft2  - -cos2itani3fj/3l , 

(Alle) 

32fj/3Q3R  ■ -tamcosAjCOS (Aj  - ft) , 

(Allf ) 

32f1/3S13V  - t^32f j/3Q3R, 

(Allg) 

32f1/3R2  - 3 2f j /3R3V  - 32fj/3V2  - 0, 

(Allh) 

32fj/3l2  * BanlS^/dl, 

(A12a) 

32fj/3x3ft  - sec2icotl3fj/3Q, 

(A12b) 

32fj/3\3R  - sec2l [-AjSlnAjSln(Aj  - ft) 

• 

+ AjCOsAjCos(Aj  - fi) ] 

(A12c) 

32f j/3\3V  - 32fj/3\3R  + tJ32fj/3x3R, 

(A12d) 
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2 ' 2 ■> 

3 fj/3ft'  - - cos*-! tan\ 3f  /3i  , 


(A12e) 


• • 

3 f j/anaR  - tarn  [AJsinA<oos(AJ  - ft)  *-  A^osAjSinCA  - ft)  J , (A12f) 


a^/anav  - 3‘fj/aftaR  + t1?2f  /»wr, 


(A12r) 


3‘f(/3R‘  - 32f  /3RDV  - a2f  /3V“  - 0. 


(A12h) 


Since  the  V derivatives  are.  In  general,  a factor  of  t » 0,007 
down  from  the  R derivatives  we  expect  poor  accuracy  for  the  final  value 
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